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ABSTRACT 

We study the statistical properties of the 1st year WMAP data on different 
scales using the spherical mexican hat wavelet transform. Consistent with the 
results of Vielva et al. (2003) we find a deviation from Gaussianity in the form 
of kurtosis of wavelet coefficients on 3 — 4° scales in the southern Galactic hemi- 
sphere. This paper extends the work of Vielva et al. as follows. We find that 
the non-Gaussian signal shows up more strongly in the form of a larger than 
expected number of cold pixels and also in the form of scale-scale correlations 
amongst wavelet coefficients. We establish the robustness of the non-Gaussian 
signal under more wide-ranging assumptions regarding the Galactic mask ap- 
plied to the data and the noise statistics. This signal is unlikely to be due to the 
usual quadratic term parametrized by the non-linearity parameter f NL . We use 
the skewness of the spherical mexican hat wavelet coefficients to constrain /jvl 
with the 1st year WMAP data. Our results constrain f^L to be 50 ± 80 at 68% 
confidence, and less than 280 at 99% confidence. 

Subject headings: Cosmology: Cosmic Microwave Background, methods: statis- 
tical 



1. Introduction 



The current cosmological model assumes Gaussian initial conditions, created by in- 
flation. This assumption regarding the nature of primordial density perturbations can be 
verified by studying the distribution of temperature fluctuations in the cosmic microwave 
background (CMB). While the simplest inflationary models predict Gaussian primordial per- 
turbations, there are other models of inflation, such as those involving multiple scalar fields, 
features in the inflaton potential or phase transitions, that could give rise to detectable non- 
Gaussianity. Hence studies of Gaussianity help distinguish between different early universe 
scenarios. Gaussianity is also a key underlying assumption of CMB data analysis wherein 
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the angular power spectrum fully specifies its statistical properties, and must be tested. 
Non-Gaussianity can also be associated with secondary anisotropies in the CMB, or with 
foreground contamination and systematic effects. 

Prior to the release of WMAP data there was no clear evidence of cosmological non- 
Gaussianity. Since the release of 1st year WMAP data, a number of tests of non-Gaussianity 
have been performed, with somewhat differing results. Each statistic is sensitive to a dif- 
ferent kind of non-Gaussianity, hence there is need for a wide variety of tests. Komatsu et 
al. (2003) use an optimized test based on the bispectrum, as well as Minkowski functionals, 
while Colley & Gott (2003) study the genus, and both groups report consistency with Gaus- 
sianity. Gaztanaga & Wagg (2003) do a 3-pt angular correlation function analysis and find 
consistency with Gaussianity as well. Chiang et al. (2003) perform a study of the phases 
of spherical harmonics and find some evidence for non-Gaussianity at high multipoles. Copi 
et al. (2003) find some evidence for low I correlations and deviation from isotropy. Park 
(2003) find a large difference between the genus amplitudes of the northern and southern 
hemispheres and a positive genus asymmetry in the southern hemisphere. Eriksen et al. 
(2003a) compute the 2 and 3-pt correlations and report a significant north-south asymme- 
try; Eriksen et al. (2003b) use Minkowski functionals and find a significant genus in the 
northern hemisphere and again indications of north-south asymmetry. Hansen et al. (2004) 
use local curvature and find non-Gaussianity/asymmetry in the data on scales of a few de- 
grees. Gurzadyan et al. (2004) find ellipticity in the temperature anisotropy features in the 
data, consistent with what was found previously in BOOMERang data. Vielva et al. (2003; 
hereafter V03) report a non-Gaussian signal in the southern hemisphere at high significance 
in the form of kurtosis on ~ 4° scales using the spherical mexican hat wavelet transform 
on WMAP data. Some of the detections of non-Gaussianity and/or asymmetry thus far 
reported in the WMAP data are at the level of 99% or greater. 

Wavelet transforms are useful tools in non-Gaussianity studies because they enable the 
signal on the sky to be studied on different scales, with simultaneous position localization, 
so that the obscuring effects of the central limit theorem, that can exist in both real and 
Fourier spaces, are reduced. With wavelets any non-Gaussian detection can be localized on 
the sky in scale and position, so that its nature and source can be better determined. Planar 
wavelets have been used in Gaussianity studies of the CMB by Pando et al. (1998), Hobson 
et al. (1999) and Mukherjee et al. (2000), while Barreiro et al. (2000) use the spherical Haar 
wavelet, and Cayon et al. (2001,2003), Martinez-Gonzalez et al. (2002), and V03 use the 
spherical mexican hat wavelet (SMHW). Wavelet methods have been compared with other 
pixel or Fourier based methods in Hobson et al. (1999), Aghanim et al. (2003), Cabella et 
al. (2004), and the performance of isotropic as well as highly anisotropic multi-scale bases 
in distinguishing between different sources of non-Gaussianity in the CMB has been studied 
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in Starck et al. (2003). 

In this paper we use the spherical mexican hat wavelet transform to probe non-Gaussianity 
in the WMAP data. We extend the work of V03 by performing new multiple tests of the 
robustness of the non-Gaussian signal. We also look at the non-Gaussianity in terms of 
an excess in the number of cold pixels, and in terms of scale-scale correlations amongst 
wavelet coefficients. Further, we place constraints on a popular form of non-Gaussianity (a 
quadratic term in the curvature perturbations parametrized by the non-linearity parameter 
Jml)- This paper is organized as follows. In §2, we present results from using the SMHW 
transform on WMAP data. Whilst confirming the results of V03, we perform new multiple 
tests of the robustness of the non-Gaussianity signal in the kurtosis spectrum (a) through 
the use of different (extended) masks, and (b) relaxing the assumption of a simplified noise 
model. We find that the signal shows up even more significantly in the form of the number of 
cold pixels (or coefficients). In §3, we examine scale-scale correlations amongst the wavelet 
coefficients. We find significant deviations from Gaussianity, a corroboration of the signal 
detected and described in §2. In §4 we obtain constraints on the non-linearity parameter 
Inl- We conclude in §5. 



2. Skewness and Kurtosis of Wavelet Coefficients 

A non-Gaussianity detection in the 1st year WMAP data was reported by V03. Applying 
the SMHW transform to the Q-V-W coadded data, and computing the skewness and kurtosis 
of the wavelet coefficients over scales ranging from about 10 arcmins to 10 degrees, they found 
that the kurtosis of wavelet coefficients on scales ~ 4° was too high at the 99.9% confidence 
level. It was found that the excess kurtosis was in the southern hemisphere, while the kurtosis 
signal in the northern hemisphere was consistent with Gaussianity. The signal was shown to 
be independent of frequency. 

It is important to determine whether the significance of this non-Gaussianity detec- 
tion is affected by systematic effects, such as the choice of mask, or simplified assumptions 
about noise. Since the spherical mexican hat wavelet transform is a sensitive probe of non- 
Gaussianity, we use it to perform an independent analysis of the 1st year WMAP data. 
The basic steps followed in the analysis are as follows: starting with the foreground cleaned 
Q-V-W coadded data, we bring the map down in resolution to Healpix n si d e = 256, apply 
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the kpO 1 mask, perform SMHW transforms 2 to obtain wavelet coefficients corresponding to 
the different scales R (also setting the monopole and dipole of the map to zero here), ap- 
ply appropriately extended versions of the mask to the wavelet coefficients of each scale to 
exclude coefficients contaminated by the mask and known point sources, and compute the 
skewness and kurtosis of the remaining unmasked coefficients. 

For our results to be directly comparable we perform the SMHW analysis for the same 
scales used by V03. For convenience and clarity, the scales Ri (i=l,2,...,15) plotted in the 
subsequent figures are listed in Table f. 



Table 1 

The scales used in the spherical mexican hat wavelet transform 



scales 


Ri 


R2 


R 3 


i?4 


R 5 


Rq 


i?7 


i?8 


Rq 


Rio 


Rn 


Rl2 


Rl3 


i?14 


Rl5 


arcmin 


14 


25 


50 


75 


100 


150 


200 


250 


300 


400 


500 


600 


750 


900 


1050 



Fig. f shows our results, with the extended mask at each scale obtained by extending 
the /cpO+sources mask such that all pixels closer than 2.5R to any of the pixels in the 
kpO+sources mask within |6| < 25° are excluded from the analysis, attempting to follow 
closely the procedure of V03. The mask was not extended around point sources outside of 
this region as the fcpO+sources mask around point sources seems to be in general extended 
enough to not cause contamination in wavelet coefficients on small scales, and on larger scales 
the effect gets averaged out. 3 Conservatively extending the region around point sources too 
would leave too few pixels on scales of interest here. (It is unlikely that the non-Gaussianity 
signal found below, on ~ 4° scales, is coming from point sources; see also Fig. 4(b) and 
related text in Sec. 2. 1). 

The mean, la, 2a and 3a confidence contours obtained from 1000 Gaussian simulations 
processed in the same way as the data are also shown in Fig. 1. The Gaussian simulations 



1 By kpO we mean the kpO mask without sources. 
2 See Appendix for further details. 

3 There are however 3 bright semi-point sources in the northern hemisphere that are seen to visibly 
contaminate wavelet coefficients on scales less than R7. These are taken care of by extending the mask out 
to 2.5i? around these sources for scales less than R^. There is also a more diffuse spot in the kpO mask 
outside the |6| < 25° region which causes visible contamination even on larger scales. The mask around this 
region is extended as well. The extended masks around these 4 regions are always retained even when the 
remaining details of the mask are varied. For reference, these regions are shown mapped in fig. 2 of Eriksen 
ct al. 2004. The actual number of pixels in the coadded map that contain emission from these 4 sources 
together is quite small (less than 40, for n si d e — 256). 
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were created in the following way. CMB realizations with the same C\ spectrum as the 
flat-A CDM cosmology with power-law primordial power spectrum that best fits the 1st year 
WMAP data (Spergel et al. 2003; Hinshaw et al. 2003) were created at Healpix resolution 
n S ide — 512. Each realization was copied and smoothed with the WMAP beam window 
functions for each of the Q, V and W radiometer channels. Independent noize realizations 
of rms do/ V N b s were then added to the maps, where the effective number of observations 
N obs varies across the sky, and a is different for each radiometer channel. The 8 maps thus 
produced were coadded weighted by N obs /al. Thereafter the same analysis procedure that 
is applied to the data map is applied to each of the Gaussian simulations. 

If we estimate the significance of the signal using the generic y 2 test, that includes 
information on all scales, 

x 2 =J2 i s w - smvjinisiRi) - siRj)}, (i) 

where S(Ri) is the skewness or kurtosis signal on scale Ri, S(Ri) is the mean obtained from 
Gaussian simulations, and the covariance matrix obtained from simulations, we find 

that by comparing the \ 2 °f the signal in the data with the distribution of x 2 's obtained 
from Gaussian realizations, the kurtosis signal in 9 of 1000 realizations have larger x 2 's than 
the data. Hence we arrive at a significance of 99% for this signal. The skewness signal is 
consistent with Gaussianity. Significances obtained using the x 2 test are tabulated in Table 
2, under mask 1. 

We find good consistency with the signal reported in V03. The kurtosis signal on 3 — 4° 
scales in the southern Galactic hemisphere is outside the 3a confidence contour. Only 2 and 
3 simulations out of 1000 lead to a stronger kurtosis signal on scales R7 and R 8 respectively, 
and only 1 of the Gaussian realizations has a larger kurtosis than the data on both the 
R 7 and R 8 scales in the southern hemisphere. Hence the signal on these particular scales 
appears significant at the 99.9% confidence level. 

Even if the data were Gaussian there is a certain probability of obtaining outlier signals 
in its kurtosis spectrum on at least two of the 15 scales considered. We take a closer look 
at this probability to better understand the significance of the detection. The number of 
Gaussian realizations that have kurtosis values that fall outside the 99% confidence region 
in the southern hemisphere in any two of the 15 scales is 28 (this number is 17 for a pos- 
itive kurtosis), indicating that the signal found above on scales R 7 and R 8 in the data is 
significant at at least the 97% confidence level. (The two scales were always consecutive 
in the simulations, but well spread out amongst all the scales.) Further these numbers are 
22 and 17 for the northern hemisphere, and in none of these cases did the signal on both 
hemispheres lie outside their 99% limits implying that a significant north-south asymmetry 
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in the kurtosis signal in any two scales was seen in 50 (37 for positive kurtosis) of 1000 
Gaussian realizations. Thus the north-south asymmetry itself appears significant at at least 
the 95% level. 

Histogram plots of the wavelet coefficients on scale R-j are shown in Fig. 2, for the all sky, 
northern and southern Galactic hemispheres. A longish tail towards negative values is seen 
in the southern hemisphere. Fig. 3 shows statistics relating to the minima, maxima (both in 
units of a on each scale) and a of the wavelet coefficients in the southern hemisphere, against 
scale, in the top panel. We use a to denote the rms dispersion of the wavelet coefficients 
on each scale, noting that la, 2a, etc., may not correspond to the same confidence levels 
as for a Gaussian distribution. The middle and bottom panels show the statistics relating 
to the number of wavelet coefficients that were larger than (mean+lcr), (mean+2<7) and 
(mean+3<r) and smaller than (mean— a), (mean-2cx) and (mean-3cr), respectively, again in 
the southern hemisphere. While the rest are seen to be quite consistent with limits obtained 
from Gaussian simulations, we see that the minima on ~ 4° scales is significant (with only 
~ 1% of the simulations showing a stronger minima on each of scales Rs and Rg), and the 
number of wavelet coefficients that are smaller than (mean-3<r) is very significant, with none 
of the simulations showing a stronger deviation on scales R$ and R?. This last estimator 
clearly gives a very strong signal. Only 3 of 1000 Gaussian realizations give a value for this 
estimator that is larger than the 99% confidence contour in any 4 scales. The signal we have 
here thus appears more significant than 99.7% as in the data the value of this estimator on 
4 scales is well out of the 99% confidence region. The x 2 test described earlier also gives a 
similar significance for this signal. The number of cold pixels in the southern hemisphere 
is too large, on scales of 3 — 4°. On scales Rq and R7 there is more than one cold spot 
contributing to this number, while on larger scales it is mainly the one cold spot located 
near (b = —57°, / = 209°) pointed out in V03. This spot is present on scales Rq and R7 
as well. The corresponding statistics for the northern hemisphere are well consistent with 
Gaussianity. 

Hence the non-Gaussianity shows up in the southern hemisphere in the form of a kurtosis 
signal and a larger than expected number of cold pixels. 

Table 2 

Significance of deviation of the skewness and kurtosis signals from Gaussianity using the \ 2 test 



masks 


mask 1 


mask 2 


mask 3 


mask 4 


ILC 


skewness 


44% 


13% 


57% 


79% 


24% 


kurtosis 


99% 


99% 


98% 


99% 


95% 
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2.1. Other Masks 

We have checked that the above results are unaltered upon using a variety of different 
masks. We now show the kurtosis signal in the southern hemisphere for a few different 
masks. The mean, la, 2a and 3a confidence contours are obtained from Gaussian simulations 
processed each time in the same way as the data. 

Fig. 4(a) show the result of using an extended mask that is made in the same way as 
for Fig. 1 above but by extending the boundaries of the whole kpO (without sources) mask 
by 2.5R, and then adding the mask around point sources back in on scales smaller than R 7 . 
The shape of the mask is then retained on all scales. 4 

In Fig. 4(b) we show the result of using extended masks that on each scale apply a 
straight |6| < (25 + 2.5R)° Galactic cut, as well as a mask around point sources for scales 
smaller than R 7 . The signal is thus unaffected by the shape of the mask. Fig. 4(c) shows the 
result of this time using a straight \b\ < (35 + 2.5/2)° Galactic cut. In going from Fig. 4(a) 
to 4(c) more of the sky is being excluded by the extended mask, and it is seen in the form 
of increased variance (this effect being larger for larger scales). But while the significance 
of the kurtosis signal seems to go down in this way, Fig. 5 shows that the number of cold 
pixels are in fact equally or more significant for this last mask. 

Fig. 4(d) shows results from the ILC map. In this case, since there is little contamination 
from the Galactic plane left in the map, we can use just the /cpO+sources mask without any 
extensions and apply it after the SMHW transform. The Gaussian simulations here were 
obtained by simulating the signal in each of the 10 radiometer channels, then smoothing to 
1° resolution, obtaining the noise weighted averaged signal for the 5 frequency channels, and 
then taking a linear combination of these with weights given in § 4 of Bennett et al. (2003). 
The signal is found in the ILC map too. 

From the above analysis it is found that the kurtosis signal is indeed independent of 
the properties of the mask. The significances derived using the \ 2 t es t f° r above masks are 
tabulated in Table 2. The masks corresponding to Figs. 4a,b,c are labelled mask 2, mask 3 
and mask 4 respectively in the table; the ILC case with the unextended /cpO+sources mask 
is labelled ILC. 



4 We have also checked that applying the kpO mask before or after performing the SMHW transform (but 
before applying the extended masks) does not affect the signal. Also, if we apply the SMHW tranform to 
the kpO mask, and identify the wavelet coefficients that are say more than 1% contaminated (by the edge of 
the mask), and include these pixels to make the extended masks, we have checked that the signal remains 
unaffected. 
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2.2. 



Noise Simulations 



Finally Fig. 6 shows the kurtosis signal for the all sky case for the same mask as used 
in Fig. 1 but this time with confidence contours obtained from 110 Gaussian simulations 
that make use of the 110 full noise simulations provided by the WMAP team for each ra- 
diometer channel. The full noise simulations were made by generating one year of simulated 
time-series data including white noise, 1/f noise, and all inter-channel correlations that are 
known to exist in the radiometers, and then taking this data through all the steps of pro- 
cessing such as flight calibration, map-making, and filtering pipeline. We found that when 
compared with the case of 110 realizations of simple white noise for the same underlying sky 
simulations the two cases give identical results, so that when plotted simultaneously they 
are indistinguishable. This indicates that our simulations are reliable, and that the simple 
white noise model is completely satisfactory according to this statistic. 



Having obtained the wavelet coefficients of the data on different scales, we compute the 
scale-scale correlations between corresponding coefficients, 



w(Ri,x) are wavelet coefficients on scale Ri, and position or pixel x in the sky. The coeffi- 
cients that contribute to the sums are the N unmasked coefficients on the larger scale. 

The top panel of Fig 7 shows scale-scale correlations between wavelet coefficients. The 
exact order of plotting is given in Table 3. Coefficients over the whole sky are used here. 
The mean, la, 2a and 3a confidence contours obtained for the scale-scale correlations from 
Gaussian simulations are also shown. Consistency with Gaussianity is indicated, amongst 
these well separated scales. 

Zooming into the scales that indicated non- Gaussianity in the previous section, the 
bottom panel of Fig. 7 shows the scale-scale correlations between coefficients of scales R e 
to R n in the southern Galactic hemisphere. The exact order of plotting is given in Table 3. 
Compared to Gaussian simulations we now find significant scale-scale correlations between 
the wavelet coefficients, presumably due to the cold spot that is common to all of these 
scales. A similar result is found for different masks. 



3. 



Scale-scale Correlations 




NY, x w{Ri,x) 2 w{Rj,x) 2 



(2) 



Table 3 
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The scales associated with the dummy indices plotted in Fig. 7 



index# 


1 


2 


3 


4 


5 


6 


7 


8 


9 


10 


11 


top panel 


Ri,R 2 


Ri,R 3 


Ri,R 6 


Ri,R 7 


ill, ilio 


Ri, R14 


R 2 ,R 3 


R 2 ,Rr, 


-R 2 ,Rt 


-R2 , Rio 


-R2 , Rl4 


index# 


12 


13 


14 


15 


16 


17 


18 


19 


20 


21 




top panel 


R 3 ,Rs 


R 3 ,Rt 


R3,Rio 


R 3 ,Rl4 


R 5 ,Rt 


R 5 ,Rio 


-Rs , R14 


R 7 ,Rio 


R7,Rl4 


Rio, -R14 




index# 


1 


2 


3 


4 


5 


6 


7 


8 


9 


10 


11 


bottom panel 


Re,R7 


R$, Rs 


Rq , Rg 


-R6: RlO 


-R6, Rll 


Rt,R& 


R 7 ,Rg 


R7, i?io 


«7,Rll 


R s ,Rg 


-Rg, Rio 


index# 


12 


13 


14 


15 
















bottom panel 


Rs,Ru 


Rg, Rio 


Rg, ilii 


Rio,Ru 

















4. Jnl Constraints 

We showed in §2 that the skewness signal obtained using a range of masks was consistent 
with Gaussianity. In this section, we use the skewness signal on the 15 scales to compute the 
limits that can be placed on the amplitude of primordial non-Gaussianity, as parametrized 
by the non-linear coupling parameter, f NL . Since the non-linear term (defined below) is not 
dominant and is quadratic, it is the skewness signal in the data, rather than any higher order 
cumulants, that is expected to be most sensitive to it. The kurtosis signal is not sensitive to 
this kind of non-Gaussianity (Cayon et al. 2003). 

/atl characterizes the amplitude of a quadratic term added to the curvature perturba- 
tions, 

= $l(x) + f NL [S£(a;) - (<£>!(*)>] , (3) 

where $l are Gaussian linear perturbations with zero mean. Thus J'nl parametrizes the 
leading order non-linear corrections to the primordial (curvature) perturbations. The mo- 
tivation to use data to place constraints on /jvl is to address how Gaussian current CMB 
data are, or how much primordial non-Gaussianity, of this particular form, is allowed by the 
data. Such analyses also help compare the sensitivity of different data sets and of different 
estimators of non-Gaussianity to this particular form of non-Gaussianity. 

Using an optimal estimator of non-Gaussianity based on the bispectrum, namely the 
cubic statistic, Komatsu et al. (2003) place limits of —58 < /jvx < 134 at 95% confidence 
using the WMAP lst-year data. They derive 95% confidence limits of f NL < 139 based on 
Minkowski functionals. Smith et al. 2004, using VSA data, obtain an upper limit of 3100 
at 95% confidence (their limit being large because their data is sensitive to small scales). 
Previously, using COBE DMR data, the bispectrum analysis of Komatsu et al. (2002) placed 
a limit of |/jvl| < 1500, and using the skewness of SMHW coefficients Cayon et al. (2002) 
placed a limit of |/jvl| < 1100, both at 68% confidence. Using MAXIMA data Santos et al. 
(2002) placed a la limit of \f NL \ < 950. 

We use the non-Gaussian simulations of Komatsu et al. (2003). How these have been 
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produced are described in detail in the Appendix of their paper. Since producing these maps 
is a computationally intensive process, we use the 300 available realizations of non-Gaussian 
sky maps at HEALPix resolution n s ^ e =256 to obtain the mean values of skewness at each 
of the different scales, for different values of Jnl- We use Gaussian simulations to estimate 
the covariance matrix of the skewness values for the different scales, and thus in turn to 
estimate the uncertainty in the measured /atl- Since Gaussian simulations can be computed 
several orders of magnitude in time faster, we can estimate the covariance matrix accurately 
using a larger number of Gaussian simulations, and the uncertainty estimated from Gaussian 
simulations is a good approximation to that estimated from non-Gaussian simulations for 
\f NL \ < 500 5 . 

In order to estimate the maximum likelihood f^L, we compare the skewness values of 
the data at the 15 scales with simulations and use the goodness of fit statistic 

X 2 = E l S W - § sim (Ri)] [S(Ri) - SsimiRj)] , (4) 

where S(R,j) is the skewness of WMAP data on the ith scale Ri, S S i m (Ri) is the mean 
value from Monte Carlo simulations, computed for different values of Jnl, and S^^^. is the 
scale-scale skewness covariance matrix from simulations. 

We have tested that the x 2 statistic accurately recovers /jvz, by using it on simulated 
maps. This is illustrated in Figure 8. The figure shows the mean x 2 distributions obtained 
from 300 simulated realizations of non-Gaussian maps with Jnl values of 50, 100 and 150. 
The simulations include noise and window functions of the WMAp 1-yr data in the same as 
described earlier. 

A plot of x 2 values obtained using data for different /tvl is shown in figure 9. 6 J'nl is 
thus estimated to be 50 ± 80 at 68% confidence, with 95% and 99% upper limits of 220 and 
280 respectively. 



5 One way to see this is that the la uncertainty indicated by the different curves in figure 8 below are 
roughly the same indicating that at current sensitivity of the data the uncertainties do not depend on the 
precise best fit value of Jnl- 

6 This x 2 plot is obtained from using a diagonal covariance matrix of the skewness values on different 
scales. 1000 Gaussian simulations may be too small to obtain convergence for the off-diagonal elements of 
the covariance matrix. However, as also discussed in Erikscn ct al. 2003b, if we are interested in obtaining 
the probability of the data given the Gaussian hypothesis, as was done in §2 of this paper, valid results 
can be obtained if we compute both the x 2 of the data and of the Gaussian realizations using a diagonal 
covariance matrix. Similarly a diagonal covariance matrix can be used here when we are interested in the 
relative change in \ 2 w hh respect to a parameter. Although we do obtain consistent limits on f NL upon 
using the full covariance matrix. 
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The limits on f^L can also be checked using the fisher discriminant function (Barreiro 
& Hobson 2001, Cayon et al. 2003). An optimal linear function of the measured variables 
(here, skewness of wavelet coefficients on 15 different scales) is 

f(x) = (/i - /ixfW^x. (5) 

Here x is the 15 element vector that contains the skewness values at the different scales 
considered here, and i(x) is the fisher discriminant function that optimally puts together 
information contained in x in the sense of maximizing the difference between the expected 
mean values of t from Gaussian and non-Gaussian models, and minimizing their dispersions. 
W = V + Vi, the sum of the covariance matrices of the test statistic in the Gaussian 
(subscript 0) and non-Gaussian (subscript 1) cases. ^ an d A^i are vectors containing the 
mean values of the test statistic in the Gaussian and non-Gaussian cases respectively. 

Thus for different Jnl, the fisher discriminant function is found for each of the Gaussian 
realizations, for each of the non-Gaussian realizations, and for the data. The probability 
that the data are drawn from one or the other hypothesis can then be estimated. Looking 
at the fraction of non-Gaussian simulations that have values larger (for positive /atl) and 
smaller (for negative /atl) we deduce similar limits on /jvx as derived above using the x 2 
test. Barreiro & Hobson (2001) found that the fisher discriminant can do better than x 2 's 
at distinguishing between Gaussian and non-Gaussian hypotheses. However since we have 
only 300 non-Gaussian simulations to obtain the fractions from, the accuracy of this method 
is not better but just comparable for our case here. Figure 10 shows the histograms of the 
fisher discriminants of 1000 Gaussian realizations, of 300 non-Gaussian realizations and of 
the data, for f NL values of 120 and 250; these values are close to the the la and 2a limits 
derived using the x 2 above. For these values of f NL , 0.74 and 0.95 of the non-Gaussian 
simulations, respectively, have larger values of the fisher discriminant than the data. We see 
that results obtained using fisher discriminants are consistent with those obtained using the 
X 2 test 7 . 

We note that if the skewness signal in the data showed deviations from Gaussianity on 
particular scales then this method involving wavelet transforms could be used to obtain scale 
dependent constraints on f^L- However the skewness signal in the data is well consistent with 
Gaussianity. The skewness spectrum of Fig 1(a) really flattens out with the use of another 
mask, such as the one used in Fig 4(a), while the kurtosis signal remains unchanged. Thus we 



7 For reference, the statistical power of the fisher discriminant test, as defined for example in Barreiro & 
Hobson (2001) is found to be 0.28, 0.66, 0.88, and 0.98 for f NL values of 120, 250, 350, and 500,respectively, 
at the 95% confidence level. These fractions are approximate as they have been obtained using only 300 
non-Gaussian realizations. 



- 12 - 



do not believe there is reason to distinguish between different scales in obtaining constraints 
on f NL here. Also, the constraints on f NL derived here could possibly be made more stringent 
if we used the above method on the Wiener filtered map of primordial perturbations as 
discussed in Komatsu, Spergel & Wandelt (2003). We will explore this in a future paper. 
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5. Conclusions 

We have analysed the first year WMAP data using a spherical mexican hat wavelet. We 
detect non-Gaussianity at ~ 99% significance, consistent with that reported by Vielva et al. 
(2003). This detection corresponds to a positive kurtosis and to the presence of a larger than 
expected number of cold pixels (wavelet coefficients) in the southern Galactic hemisphere on 
scales 3 — 5°. 

We have tested for changes in the significance of the signal with the type of mask used. 
The signal is found to be robust, and is found in the ILC map as well. We have also compared 
confidence contours obtained for the kurtosis spectra using the full noise simulation maps 
provided by the WMAP team, containing 1/ / noise and other effects from data processing, 
to those obtained from using simulations that contain just white noise. We find very good 
agreement. 

We have also applied another test statistic, the scale-scale correlations between wavelet 
coefficients. Significant scale-scale correlations are seen amongst the coefficients over the 
range of scales that indicate the above non-Gaussianity. 

We then use the skewness statistic on the different scales to place constraints on the 
non-linear coupling parameter f NL , the motivation being to see how much non-Gaussianity 
of this particular form is allowed by current data. It is also a way to compare the sensitivity 
of different test statistics to this parameter. Constraints obtained are closely consistent with 
those obtained by Komatsu et al. (2003) using the cubic statistic and Minkowski functionals 
on the same data. The constraints on /jvl derived here could possibly be made more stringent 
if we used spherical wavelets on the Wiener filtered map of primordial perturbations as 
discussed in Komatsu, Spergel & Wandelt (2003). We will explore this in a future paper. The 
kurtosis statistic is not sensitive to this form of non-Gaussianity. We will present constraints 
on other forms of non-Gaussianity implied by the kurtosis statistic of the WMAP data 
elsewhere. 
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Fig. 1. — Skewness and kurtosis spectra obtained using a SMHW analysis of the WMAP 
coadded data (stars), together with the mean values, and la (solid), 2cr (dashed) and 3cr 
(dotted) confidence contours obtained from 1000 Gaussian simulations. The top panel shows 
these for all sky, the middle for the northern Galactic hemisphere, and the bottom panel 
shows these for the southern Galactic hemisphere. Similar results were presented in V03. 
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Fig. 3. — Statistics relating to the minima and maxima, in units of a, and a itself (top panel), 
and the number of wavelet coefficients that are larger than (mean+cx), (mean+2<r) and 
(mean+3<r) (middle panel), and smaller than (mean— a), (mean-2a) and (mean-3cr) (bottom 
panel), are shown, all for the southern Galactic hemisphere. Comparing the data (coadded 
WMAP; stars) with the confidence limits obtained from Gaussian simulations (mean and 
la: solid, 2a: dashed, and 3a: dotted), the data show an excess of cold coefficients. 
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Fig. 4. — Kurtosis spectra of wavelet coefficients in the southern Galactic hemisphere for 
different masks, as described in text. 




Fig. 5. — Statistics relating to the minima, shown in units of a, and to the number of 
wavelet coefficients that are smaller than (mean-3cr), in the WMAP coadded data (stars) in 
the southern Galactic hemisphere for a mask that is a straight |fe| < (35 + 2.5i?)° Galactic cut 
(see text) are shown. Also shown are the mean and la (solid), 2a (dashed) and 3a (dotted) 
confidence contours obtained from Gaussian realizations. 
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7. Appendix 



The SMHW transform has been used previously by Cayon et al. (2001,2003); Martinez 
Gonzalez et al. (2002), and V03. It is a continuous and symmetric wavelet which in the 
small angle limit corresponds to the euclidean mexican hat wavelet. The wavelet is given by 



! + (§ 



1 2 



y^2 



-y 2 /2R 2 



/2ttN(R) 

where R is the scale and N(R) is a normalization constant given by 

N(R) = R(1 + ^ + 



(6) 



(7) 



y = 2tan^, for polar angle 6. ip is thus isotropic so that when centered on the north pole 
it is independent of the polar coordinate 0. Thus moving any point to the north pole, the 
wavelet coefficients 

(8) 

can be obtained via a convolution of the sky map T(x : 0, 0) with the wavelet function. This 
is easily performed in spherical harmonic space. 



w(R,x' : (9', (/>'))= / dxT(x)V R (\x? -x\) 
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Fig. 6. — This plot is to be compared with Figure 1(b). The confidence contours here have 
been obtained from 110 simulations using the 110 full noise simulations for each radiometer 
channel made available by the WMAP team. 
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Fig. 7. — The top panel shows scale-scale correlations amongst wavelet coefficients between 
scales that are well separated and span the whole range studied here. The bottom panel 
looks at scale-scale correlations amongst wavelet coefficients in the southern Galactic hemi- 
sphere between scales Rq to Rn, the scales that show excess kurtosis and an excess of cold 
coefficients. See Table 3 for the exact order of plotting. The correlations obtained from 
the WMAP coadded data are shown as stars. Also shown are the mean and la (solid), 2a 
(dashed) and 3a (dotted) confidence contours obtained from Gaussian realizations. 
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Fig. 8. — Mean x 2 distributions obtained from 300 simulated realizations of non-Gaussian 
maps with f NL values of 50 (dashed), 100 (solid) and 150 (dot-dashed), illustrating that 
the statistic recovers the input value of /nl- The simulations include as usual the noise 
properties and window functions of the WMAP 1-yr data. 




Fig. 9. — A plot of x 2 against f^L obtained using WMAP data. f NL is thus estimated to be 
50 ± 80 at 68% confidence, and the 95% and 99% upper limits are 220 and 280 respectively. 




Fig. 10. — Histograms of fisher discriminant values from 1000 Gaussian simulations (solid), 
from 300 non-Gaussian simulations (dashed), and the fisher discriminant of the data (vertical 
line), for f NL values of 120 and 250. These values are near the la and 2a limits on f NL . 



